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Asymptotic preserving schemes on distorted meshes for 

Friedrichs systems with stiff relaxation: application to 

angular models in linear transport. 

Christophe Buet* Bruno Despres^& Emmanuel Franck* 
April 11, 2013 

Abstract 

In this paper we propose an asymptotic preserving scheme for a family of Friedrichs systems 
on unstructured meshes based on a decomposition between the hyperbolic heat equation and a 
linear hyperbolic which not involved in the diffusive regime. For the hyperbolic heat equation 
we use asymptotic preserving schemes recently designed in [FHSNII]-[BDFI1]. To discretize 
the second part we use classical Rusanov or upwind schemes. To finish we apply this method 



^^ for the discretization of the Pn and Sn models which are widely used in transport codes. 
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1 Introduction 

We study the finite volume discretization of general linear hyperbolic systems with stiff source 
terms depending of a relaxation parameter e. which admit an asymptotic diffusion limit. This 
type of system occurs in many physical applications (transport of particles, damped waves, elec- 
tromagnetism, linearized gas dynamic, plasma physics) or in biology, and poses some numerical 
difficulties. The classical Godunov-type discretizations (upwind, Rusanov or HLL schemes) are not 
efficient because the time and spatial steps are constrained by the relaxation parameter e [BDF11], 
[JL96], [Jin99]. To treat this problem S. Jin, C. D. Levemore [JL96]-[Jinll] using the ideas of A. 
Y. Leroux [GL96], introduced the notion of asymptotic preserving schemes (AP schemes) which 
eliminate these constraints. To illustrate the advantage of asymptotic preserving discretizations, 
we propose a simple numerical example. We solve the hyperbolic heat equation 



d t p+ -d x u = 0, 

£ 

d t u H — d x p H — u = 0, 

£ £ 



(1) 



with two schemes: the upwind scheme and the asymptotic preserving scheme [GT01]. This model 
is approached when e is small by the following diffusion equation 



dtp- d x ( -d x p 



0. 



The initial data is given by p(x, t = 0) = G(x) with G(x) a Gaussian function and u(x, t = 0) = 0. 
The parameters are given by a — 1 and e = 0.001. The time discretization is explicit and the time 
step is the half of the stability limit time step. The convergence errors are computed using the 
exact diffusion solution. The results proposed in table (1) and on figure (1) show that asymptotic 
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Figure 1: On the left: numerical solution of the Gosse-Toscani scheme for 50 and 500 cells, on the 
right: numerical solution of the upwind scheme 500, 1000 and 10000 cells 



preserving scheme is more precise and cheaper in CPU time than the classical upwind scheme. 
These remarks may justify to use asymptotic preserving for this type of problem. 
In ID, many AP schemes have been designed: a non exhaustive list is S. Jin, C. D. Levermore [JL96] 
or L. Gosse, G. Toscani [GT01] for the hyperbolic heat equation, M. Lemou, L. Mieussens, N. Crou- 
seilles [LM07]-[CL11]-[CL11] for some kinetic equations, C. Hauck, R. G. McClarren [HLMclO] for 
the P/v equations, L. Gosse [Gossll], C. Buet and co-workers [BCLM02] or S. Jin and C. D. Lev- 
ermore [JL91] for S N equations and C. Berthon, R. Turpault [BCT08]-[B010]-[BLeFTll]-[BT10] 



Schemes 


L 1 error 


L 2 error 


CPU time 


AP scheme, 50 cells 


0.0065 


0.0110 


0m0.054s 


AP scheme, 500 cells 


0.0001 


0.00018 


0ml5.22s 


upwind scheme, 500 cells 


0.445 


0.647 


0m24.317s 


upwind scheme, 1000 cells 


0.279 


0.113 


2m9.530s 


upwind scheme, 10000 cells 


0.0366 


0.059 


1485m4.26s 



Table 1: Table with numerical error and CPU time associated to the upwind and Gosse-Toscani 
schemes. 



for generic systems and a non linear radiative transfer model. 

For some applications (ICF simulations [DW94]) we are interested in, the stiff hyperbolic systems 
are coupled with Lagrangian hydrodynamic codes which generate very distorted meshes. Conse- 
quently it is important to design cell-centered asymptotic preserving schemes for the Friedrichs 
systems with a valid asymptotic diffusion limit on unstructured meshes. Currently these types of 
schemes based on the nodal scheme [BDFll]-[BDFCras]-[BDFproc] or the MPFA scheme [FHSN11]- 
[BM06]-[AE06] have been only designed for the hyperbolic heat equation and a non linear system 
used in radiative transfer. 

The purpose of this paper is to extend Godunov-type asymptotic preserving schemes for the 
Friedrichs systems on unstructured meshes. Firstly we introduce the Friedrichs systems and give 
a formal proof of the existence of the diffusion limit. In the second part we define a numerical 
strategy based on a decomposition between a "diffusive" part similar to the hyperbolic heat equa- 
tion and a "non diffusive" part which is negligible in the diffusion regime. This decomposition, 
close to the micro-macro decomposition [LM07] allows to design a very simple method to discretize 
stiff hyperbolic systems. Indeed, using an asymptotic preserving scheme for the "diffusive" part 
(nodal asymptotic preserving for example [BDF11]) and a classical hyperbolic scheme for the "non 
diffusive" part we obtain an asymptotic preserving discretization for the complete system. After 
this, we show how angular discretizations such as P/v and Sn models fall within this framework. 
This could be applied to other angular discretizations like those based on wavelet expansion for 
instance. To finish we propose some considerations on temporal discretizations and numerical 
results for P/v and the Sn systems. 



2 Friedrichs systems 

2.1 Definition 

In this section we introduce linear Friedrichs systems with stiff source terms and their diffusion 
limit. We work in dimension two, D C M. is a polygonal domain. 

Definition 2.1. The sub-class of Friedrichs systems that we consider are defined by: 

d t \J + -A 1 d x \J + -A 2 d y V = -^RV, (2) 

with U : D x R+ — > W l , A\, A 2 , R are constant, symmetric and real square matrices. We assume 
moreover that the matrix R is non negative, i.e. (Rx,x) > for all x £ M™ and Keri? ^ 0. 

Non invertibility of the matrix R is important to obtain a non trivial asymptotic regime. The 
parameter a is, in general, a positive and a lower bounded function but for the theoretical analysis 
we assume that a is constant and positive. The relaxation parameter is e £ ]0, 1]. Since the 
matrices A±, A 2 are symmetric, the system is hyperbolic. Indeed the matrix A±n x + A 2 n y is 
symmetric for all n = (n x ,n y ) £ S . 

We define the functional spaces: 



L 2 (D) = {V,\\V\\ LHD) 



(U, U) dxdy < oo 



{a+b<p | 

UeL 2 (D),\\lJ\\ HHD) = Y, ||flx-, tf »U|| ia <ooL 

Foremost, we recall a classical result of stability for such systems. 

Proposition 2.1. If D — [0, l] 2 wii/i periodic boundary conditions, systems (2) are stable in 
H p (D). 



Proof. We begin by proving the L 2 stability: 

^4||U||| 2( m= [ (V,d t V)dxdy = -- f (AAV + A 2 d v V + -RU,V)dxdy. 

I at J D e J D e 

Since the matrix A\ is constant and symmetric real, the first integral of the right hand side writes 

/ (;4i0 x U, U)dxdy = I f d x (U, A x \J)dxdy, (3) 

JD * JD 

and since we consider periodic boundary conditions, J D (Aid x TJ ', U)dxdy — 0. In the same way we 
show that the second integral of the right hand side is null. Thus: 



U||i a(D) = - / ^(RU,U)dxdy<0. 



l d 

2dt j ,, 



since the matrix R is non negative. The L 2 (D) norm of the solution decreases, thus the system is 
L 2 -stable. We can check that V = d x ad y b\J is also a solution of the system (2) which gives the 
HP(D) stability. D 

2.2 Diffusion limit of the Friedrichs systems 

In this section we propose a formal existence result for the asymptotic diffusion limit . We introduce 
a structure assumption. 

Assumption (Hi): Let (Ei,...E„) be the eigenvectors of R and let (Ei,...E p ) be the basis of 
the kernel of R. The vectors are orthonormal. We assume that, there are two particular linearly 
independent eigenvectors E^ , Ei 2 associated to eigenvalues Xi 1 > 0, Xi 2 > with the structure 
assumption 

A 1 E l =7 l 1 E n V»G{l...p}, ,., 

A 2 E l=7l 2 E l2 VtG{l...p}. [V 

In other sections we will show that the simplified models as Pjv or Sn in linear transport theory 
satisfy the previous structure assumption. For the hyperbolic heat equation extended to 2D p = 1, 
AiEi = E2 and A2E1 = E3 where Ei is the eigenvector associated to the eigenvalue and E2, E3 
are the eigenvectors associated at the eigenvalue 1 of the matrix R. 

Proposition 2.2. If the assumption (Hi) is satisfied, the system (2) admits the formal diffusion 
limit 

dtV ~ ^—Kxd xx V - -r^—K 2 d yy Y = 0, (5) 

A^cr Aj 2 <7 

with V = ((U,Ei),....,(U,E p )) e RP, ifi=7 1 «7 1 eR' , x W, K 2 = 7 2 ® 7 2 € R p x R p non 
negatives symmetric matrices where the vectors j k — (7* , ...,7p) are defined in (4)- 

Proof. Using a Hilbert expansion U = U + eUi + e 2 U 2 + o(e 2 ) in (2), we obtain the hierarchy of 
equations: 

\ : i?U = 0, (6) 

£ z 

\ : AAUo + A 2 d y \J Q = -aRUi, (7) 



£ 



\ : d t U Q + A 1 d x \] 1 + A 2 d y V 1 = -<jRU 2 . (8) 



Equation (G) shows that Uo £ Ker R. Therefore, 

U = ^(U,E J )E i . (9) 

Equation (7) implies the existence of Ui up to an element of the kernel under the following 
compatibility condition 

Ai<9 x U + A 2 d y V £ (Keri?)- 1 . 

The assumption (Hi) and the equation (9) show that Uo is such that 

A t d x \J = p^Uo.E^] E 4l , A 2 d y V = p.fUo,^ E, 2 . 

Using the definition of the eigenvectors and the linearity, we obtain the relation 

D / iAUo , MdyV \ . „ TT . „ TT 



which gives the expression of U 



1 (A x d x V Q , A 2 d y U 



"' \ Ajj Xi 2 



U 1= -ipp^ + ^m +z, ,zeKcri? (10) 



Using the relation (8), we show the existence of U2 up to an element of the kernel under the 
following compatibility condition 

d t U + A 1 d x \3 1 + A 2 d y \J 1 £ (Keri?) x . 

Since Ker(-R) = Vcct(Ei, ....,E p ), the compatibility condition can be written as 

ft(Uo,E i ) + 0«(i4iUi ) E 4 )+fy(4,Ui,E i ) = O ie{l..p}. (11) 

Now we plug the relation (10) in (11) to obtain the equations 



ft(Uo,B5i) - ^-5^ (AiUoAiE,) - ^—d yy (A 2 U ,A 2 E i 
AjjCT Aj,<7 

1 „ , . -. . _ , 1 



-d xy (A 1 V ,A 2 -E i ) - -^-d yx {A 2 V 0l A{E ll )+N l = 0, for t G {l..p}, 
where 



(12) 



A i2 CT A^er 



N % = d x (A 1 z,'E i ) + d v {A 2 z,E i ) = d x (z,A 1 E i ) + d y (z,A 2 E i ). (13) 

The assumption (Hi) and the orthogonality of the eigenvectors show that the terms A\Ei, A 2 E; b 
are orthogonal to z. Consequently the terms iV$ (13) are equal to zero. Now we study the cross 
terms (AiUo,^4 2 Ej) and (A 2 U , A]Ej). One has 



;j 



(AiUcAaEi) = ( ^1 ( ^(U ,E,)E, ,A 2 E< = ^ 7 )(U ,E i ) ] E il)T fE i2 ] = 0, (14) 



v \ \ I I p 



(4,Uo,AiEi) = U 2 ^(Uo,E,)E, .AiEj = ^ 7 |(U ,E,) E^E^ =0. (15) 

The cross terms vanish. For the other terms we obtain 

(A k \J,A k E % ) = ^(Uo.EjJ^E^AfcEi) = ^(UcEjfr^tf. (16) 

So the equations (12) with Uo = U are equivalent to the equations (2) with K\ — 7 1 <E> 7 1 and 
^t" 2 = 7 2 <8> 7 2 - These matrices are symmetric by definition. Moreover 

(X, K k X) = ( 7 fe , X) 2 > V X £ R p , 

therefore the matrices K k are non negatives. □ 



Remark 2.3. • Since the matrices K\ and K 2 are non negatives, the system (5) is dissipative. 

• The size of the diffusion equation (5) is equal at the multiplicity of the eigenvalue of the 
matrix R. 

• The hypothesis (Hi) is sufficient but not necessary. The assumption AfcEj £ (Kerii) for 
i € {1-vP} is also possible. 

• // Xi t = Xi 2 the diffusion equation is isotropic. 

• In 3D the proof uses the same principle. 

3 Discretization strategy 

In this section, we propose a strategy to design finite volume schemes valid for Friedrichs system 
on unstructured meshes. The method is only valid for the Friedrichs systems which have a scalar 
diffusion limit (dimKeri? =1)- The idea is to split the Friedrichs system between a "diffusive" 
part similar to the hyperbolic heat equation and a "non diffusive" part which is negligible in the 
diffusive regime. This method is in the principle very close to the micro-macro decomposition used 
in [LM07]-[CL11]. 

3.1 Principle of the "diffusive - non diffusive" decomposition 

The "diffusive - non diffusive" decomposition uses the particular structure of some Friedrichs 
systems described by the following assumption. Assumption (H 2 ): Assume that 

dim(Keri?) = 1, consequently Ker_R = Vect(Ei), 
(H 2 ) ^ A2 = A3 = A (we study isotropic diffusion limit), (17) 

AiEj = aE 2 , A 2 B 1 = aE 3 , 

with Xi the eigenvalues , by convention Ai = 0, and E^ the eigenvectors of R. 

Since R is symmetric the matrix can be written on the following form R — QDQ 1 with D di- 
agonal matrix and Q an orthogonal matrix where the column are the eigenvectors of R. Since R 
is non negative, the coefficients of D are non negative. We define V = Q'U, 



d f y 



e 



X - e [A x d x Y + A 2 d y V^j = -^DV, (18) 



with A\ = Q f AiQ, A' 2 — Q l A 2 Q symmetric matrices. 

Lemma 3.1. Under the assumption (H 2 ), the matrices A' x , A' 2 have the following block structure 



Ci „ ( C, 



A ' a _ ( C\ B 1 ) A ' 2 ~ \ C\ B. 



where B\ and B 2 are (n— 1) x (n— 1) symmetric matrices, and C\ and C 2 are I x (n— 1) matrices 
whose elements are defined by C^, j = aSkj for k = 1, 2 and j = 1, ..., n — 1 and dkj stands for the 
Kronecker product. 

Proof. Let us consider the matrixA^ . Using the definition of the matrix Q we have 

A' 1>ij = (E i ,A 1 E j ). 

For the first line we have then, remembering that E\, ..., E n is an orthonormal basis, 

A' hlj = (Ei.^Ej) = (AiEi.Ej) = (Ei.AiEj) = (aE 2) E 3 -) = aS 2j . 

Since we are dealing with symmetric matrices, the results follows for the first column. By the same 
way we obtain the desired result for the matrix A 2 . □ 




Therefore we can rewrite the system (18) as 

d t V + 1 (P hx d x V + P h yd v V) + - e (X[d x \ + A 2 d y \) = ~^DV, (19) 

where the matrices P\ tX , Pi.y are defined, as block matrices, by 

Pi,, = { Q 1 J Pi,v 
where Q\ and Q 2 are the 3x3 matrices 

a 

(J.I. : : J a 



and A 1 = A 1 — P\ iX , A 2 = A 2 — P\ lV are symmetric matrices where the first line and column of 
the matrices A 1 , A 2 are equal to zero. Next we decompose the model (19) between two systems. 
The first part of the system is very close to the hyperbolic heat equation 

d t V + X - (Pi, B a x V + Pi, y d y V) = - J-d'v, (20) 

with, for the diagonal matrix D' , D n = 0, D 22 = D 33 = X 2 and D Li = i > 4. The second system 
is given by 

dtV + X £ (A^d x V + A 2 d y v) = -^D"V, (21) 

with, for the diagonal matrix D", D lx — D 22 — D 33 — and D u = Xi i > 4. This decomposition 
is a little bit different from the micro-macro decomposition. Indeed when we use the micro-macro 
decomposition for the linear kinetic equation (some Friedrichs systems can be interpreted as angu- 
lar discretization to the linear kinetic equation) we split the isotropic part homogeneous to O(l) 
and the residual homogeneous to 0(e). When we apply an asymptotic analysis to our decomposi- 
tion we remark that we split the equations associated to the unknowns homogeneous to 0(1) and 
0(e) (first system) which gives the diffusion limit and the equations associated to the unknowns 
homogeneous to 0(e 2 ) (second system) which are negligible in the diffusion regime. 

Principle of discretization: 

The proposed numerical method consists to use an asymptotic preserving scheme for the "diffu- 
sive " part (20) and a classical hyperbolic scheme for the "non diffusive " part (21). In the following 
section we will introduce the different numerical schemes for the two parts of the decomposition. 

3.2 Asymptotic preserving scheme for hyperbolic heat equation 

The discretization of the "diffusive" part (20) is based on a specific asymptotic preserving scheme 
that we recall now for the hyperbolic heat equation 

d t p A — divu = 0, 

(22) 
„ (x a A 
o t u+ -\7p = ^-u, 

where p € K and uel 2 . In [BDF1 1] we have observed that the classical extension of Godunov-type 
asymptotic preserving schemes (Jin-Levermorc scheme [JL96] or Gosse-Toscani [GT01] scheme) in 
2D are convergent only on regular meshes which satisfy the Delaunay condition [EGH00]. Indeed 
the numerical viscosity of the hyperbolic scheme gives a non consistent limit diffusion scheme 
(two point flux approximation (TPFA) scheme [BDF11]-[EGH00]) on unstructured meshes. To 
solve this problem two methods have been proposed. In [BDF11] the extensions of Jin-Levermore 
scheme and Gosse-Toscani scheme have been designed using the nodal finite volumes formulation 
(the fluxes are localized at the nodes) [KD10]-[CDDL09] because the numerical viscosity of this 



scheme has a better structure. Another method is introduced in [FHSN11] based on the convergent 
diffusion scheme MPFA (Multipoint Flux Approximation) [AE06]. 

Let us consider an unstructured mesh in dimension two. The mesh is defined by a finite number 
of vertices x r and cells f2j . We denote Xj a point arbitrarily chosen inside Qj . For simplicity we 
will call this point the center of the cell. By convention the vertices are listed counter-clockwise 
x r _i, x r ,'K r+ i with coordinates x r = (x r ,y r ). The length lj r and the normal n jr associated to the 
node r et the cell j are defined by 



jr — ^l x r+l — x r-l| an d li-jr 



21 



J i' 



-Vr-1 +Vr+l 



(23) 



The convention is that the norm of a vector x £ M 2 is denoted as |x|. The scalar product of two 
vectors is (x, y). The JL-(b) nodal- AP scheme (2-D extension of the Gosse-Toscani scheme) writes, 




Figure 2: Notation for node formulation. The corner length lj r and the corner normal x\j r are 
defined in equation (23). Notice that (j>n, r is equal to the orthogonal vector to the half of the 
vector that starts at x r _i and finish at x r+ i. The center of the cell is an arbitrary point inside 
the cell. 



see [BDF11], 



Qj | d t pj H — yj(^>Af r u r , n jr ) = 0, 

r 

% | d t Uj + - ^2 a jr M r (uj - u r ) = — I y^ a jr (T d - M r ) J u,-, 

r \ r / 



with the fluxes 

^2 «ir ) U r = Yl tjrPjn-jr + ^ S ^ U J ' M r = I J2 



Jr 



(T\ 



as 



Yl%A \Ya jr 



3/3 3 \ 3 3 

and the tensors 

iL'tI 7" — v i f 1A 7 t^ V_y \ -"-7" "■ 7 / ; £-t q y — t 7 7 1 A A 7 7" v-y AA«f j - * « 

This scheme admits the following limit diffusion scheme on coarse grids 

a 2 



\Qj\d t Pj(t) H Y] l jr (u r , Ujr) = 0, 

r 



(24) 



(25) 



(26) 



We recall some properties of this scheme: 



• The scheme (24) is stable for the L 2 norm [BDF11]. 

• The matrix A r = y\ Zj r iij r (8)(x T .— Xj) is positive and coercive under non restrictive conditions 
on the meshes. 

• In [BDF11] we prove that the limit diffusion scheme is convergent if the matrix A r is coercive. 

• If we implicit the source term of (24), we observe numerically that the stability CFL condition 
is independent of s. 

• These schemes exhibit spurious mods [BDF11]-[CDDL09] which degrade the quality of the 
numerical solution. A solution, based on geometrical corrections, to treat this problem is 
proposed in [BDF11]. 

• Numerical tests show convergence in all cases. 

3.3 Numerical schemes for the hyperbolic "non-diffusive" part 

We discretize the "non-diffusive" part (21) using a classical hyperbolic scheme. First of all, we 
recall two different schemes, the upwind scheme and the Rusanov scheme. The upwind scheme in 
dimension two has been studied in [Cou06] . Consider 

d t V + M t d x V + M 2 3 V \J = 

U(t = 0) = U o , ( ' 

with two arbitrary real symmetric matrices Mi and Mi. 

For a cell of index j, iijk denotes the outward normal associated to the interface dfijk between 
the cell j and one of its neighbors of index k, Gjk = M\n x , k + Min v - k = —Gkj and Ijk = |9f2^|. 

Definition 3.1. The space discretization of the upwind scheme is 

l^lftU^+^fcU^-O, (28) 

k 

with the fluxes 

Ujk = (G 3 k) + V 3 + (G jk )-U k , (29) 

where the matrices (Gjk) + '~ are the non positive and the non negative parts of the matrix defined 
by G, k ' = P~ l D + >~ P with D +, ~ the matrices of the positive and negative eigenvalues of Gjk and 
P is the orthogonal matrix such that PGjkP -1 is diagonal. 

The computation cost associated to the upwind scheme can be important for large linear system. 
Therefore we propose another choice: the Rusanov scheme. This scheme use only an estimation of 
the maximal eigenvalue. 

Definition 3.2. The Rusanov scheme is defined by 

| % | d t V j +J2ljkV jk = 0, (30) 

k 

with the numerical fluxes given by 

Uj + Ufe Uf - Ufc 
Ujfc = Gjk h Sjk , (31) 

and the local speed Sjk is such that Sjk > maxi(A* fc ) and A* fc are the eigenvalues of Gjk- 



3.4 Structure of the algorithm 

We are now ready to recapitulate the explicit version of the proposed "diffusive-non diffusive" (19) 
decomposition . 

Algorithm 3.2. 

• Preparation 

— Step 1: We diagonalize in the basis of R the system 

d t V+ 1 A 1 d x U+ 1 A 2 d y U = -^RU. (32) 

e e e 2 

— Step 2: We decompose the diagonal system 

d t V + -A[d x V + -A'd v V = -^DV, (33) 

e e e z 

with V = Q*U, A[ = Q l A x Q et A' 2 = Q l A 2 Q. We obtain 

d t V + l - {P ltX d x V + Pi, y d y V) + X - (X'g> x V + 4'flwV) = ~^DV. (34) 

— Step 3: The system homogeneous to the hyperbolic heat equation is 

d t Y + J (P x , x d x V + Pi, v d v V) = -^DV. (35) 

Using an asymptotic preserving discretization such as the JL-(b) (24)-(25) scheme or 
the P\-MPFA scheme [FHSNllj, we define a matrix P^ such that 

yn+l = V „ + At p hV n (36) 

is an explicit discretization of (35). 

— Step 4- The second system is 

d t V + * (X[d x \ + A^dyW) = -^D"V. (37) 

Using a standard finite volume scheme such as Rusanov (30)-(31) or upwind (28)-(29), 
we define a matrix A^ such that 

yn+l = v „ + AtAhV n (38) 

is an explicit discretization of (37). 

• Loop in time 

— Step 1: V™ = Q*U£ 

— Step 2: We apply the explicit scheme 

Vl+ l =Vl + M{P h + A h )Vl (39) 

— Step 3: U£ +1 = QV]l +1 

Remark 3.3. In this work we diagonalize the system to obtain the primitive variables at each time 
step. However it is not necessary, we can diagonalize the system only at the first time step. 

Remark 3.4. Usual boundary conditions are easy to incorporate in (39) with standard technics 
such as the ghost cells method. 
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4 Time discretizations 

Now we quote some remarks about the time discretization. The stability condition of the time 
scheme associated to the "diffusive - non diffusive" decomposition is given by the stability conditions 
of each part of the decomposition. Some asymptotic preserving schemes used for the "diffusive" 
part and classical schemes used for the "non diffusive" part admit CFL conditions dependent of e. 
To treat this problem, we can use a fully implicit scheme or design semi-implicit scheme stable on 
the CFL condition independent of e [BDFff]. Firstly we will study the implicit discretization of 
the "diffusive - non diffusive" decomposition. 

4.1 Implicit discretization 

We study the L? stability of the implicit of the algorithm (39). The standard L? norm is ||V||? 2 = 
Y^j |^il(Vj,Vj) and the scalar product is (U, V) = ^ |fij|(Uj, Vj). Let us assume for simplicity 
that periodic boundary conditions are used so that (X, -P/jX) < for all X <E R nx ™<= and n c is the 
number of cells: this is proved in [BDF11] for the JL-(b) scheme. Moreover one has (X, AhX.) < 
for standard upwind discretizations. 

Proposition 4.1. The implicit scheme 

MV n+l = M yn + At p h yn+1 + AtAh yn+l ( 4 q) 

is stable in the norm L 2 (D). 

Proof. By multiplying (40) by V^ +1 we obtain 

(MV£ +1 ,V£ +1 ) = (MV£, V£ +1 ) + At(P h V?+\Vl +1 ) + At(A h V™ + \V% +1 ). 



To conclude, using the inequalities (V™ +1 , PhV^ 1 ) < 0, (V£ +1 , i ft V^ +1 ) < and the Cauchy- 
Schwartz inequality we obtain 

||V" +1 || L2 < ||V"|| L2 . 

□ 



4.2 Semi-implicit schemes 

We design semi-implicit schemes modifying the "diffusive - non diffusive" decomposition to obtain 
a restrictive-less CFL. We study the scheme for the "diffusive" part. In fD the JL-(b) scheme 
(24) which is equivalent to the Gosse-Toscani scheme is stable for the L°° norm under the CFL 
condition [BDFff]: 

A (M Ma\ , 

with M — „ y , . The previous CFL condition is equivalent to 



At (1) < 1. (41) 

If we use a local-implicit discretization for the source term we obtain 

M (j^) < 1. (42) 

A reasonable CFL condition is the sum of the classical hyperbolic CFL condition and the parabolic 
CFL condition. These remarks show that we can obtain a stability condition independent of e for 
the "diffusive" part using the semi-implicit JL-(b) scheme (extension in 2D of the Gosse-Toscani 
scheme). However, for the "non diffusive" part the CFL condition of semi-implicit scheme is close 
to (41) in ID. Therefore we propose to multiply the Rusanov or upwind fluxes by an adapted 
factor M in the "non diffusive" part and use a local implicit discretization of the source term. This 

If 



strategy allows to obtain CFL condition close to (42) for the complete system. 

The factor M depends on hyperbolic system studied and the scheme used. For the Rusanov scheme 

where the Rusanov velocity is Sjk and — the diffusion coefficient, the factor M is defined by 



M 



jk 



2SjkS 
2SjkS + c <jjkh' 



with h a quantity homogeneous to the characteristic length of the mesh. For example we can use 

h = d(xj, Xfc) with Xj, X& the center of the cells. 

For the upwind scheme with a velocity Xjk and the same diffusion coefficient Mjk is defined by 

Mj k = 3 — . 

2Ajfe£ + c a jk h 

5 Applications to the P/v models 

The transport of some type of particles is described by the following transport equation with 
scattering term (for example: radiative transfer equation, neutron transport linear equation) 

9 t /(x, n, t) + -n.v/(x, n, t) = ^ / (7(x, si , t) - /(x, n, tj) dsi' . (43) 

£ e 1 J S 2 V / 

The P/v systems are obtained expanding the equation (43) on the spherical harmonics basis. By 
construction, the P/v approximation is a Friedrichs system. But simplifications for 2D flows leads 
to nonsymmetric systems. The 2D form of the P/v equations (see [Bru02]-[BH05]-[Bru05]) is 



c <kh + -rfx \-^i-\ h-\ + v i+i h+i + h i-i Vi ~*i+i h+i ) (44) 



1 +d z {AT^ir-x + sr+jr+i) ^r = o 

for m t^ and 



-9Ji + 2 d * ( E i-ih-i + Fi+ili+i) / 45 \ 

+d z {Atjli + B° l+1 lf +x ) + o (I°8 l0 If) = 0, 



for m = 0. The coefficients are defined by 



, ,. ... _ (l-m)(l + m) 

h -\l „,,u„,,, A -V(2J + 1)(2*-1)' (46) 



l(l-m+ 1)(Z + m + 1) 


Y (2Z + 3)(2Z + 1) 


/(Z + m+l)(Z + m + 2) 


(2Z + 3)(2Z + 1) 


/(Z-m + l)(Z-m + 2) 



_ ) V . , ... , ^y. , ... , -, _ / (Z-m)(i-m-l) 

Q -\/ W ,,v„,, A "V (2/ + l)(2/-l) ' (47) 



r= / v - i -/v - ; -_1 F m = (l + m)(l + m-l) 
1 V (2Z + 3)(2Z + 1) ' Y (2Z + 1)(2Z-1) ' l ; 

with A™ j = P™, C, m = P/++ 1 and P™ = E^+ l . The system formed by the equations (44)-(45) 
is not symmetric. However, by an elementary change of unknowns we obtain a symmetric system. 
If we note P/™ the unknowns defined by: 

• 1 l — 1 l ■ 

• ip = -Vzi™, 

then the P/v system associated to I™ is symmetric. 
The P/v systems satisfy the following properties 

• R is a diagonal matrix. is an eigenvalue with the multiplicity 1 and 1 is an eigenvalue with 
the multiplicity n - 1 [Bru05]-[Bru02]. 
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• The eigenvalues of the system are included in ] — 1, 1[ [Bru05]-[Bru02]. 

• The hypothesis (H 2 ) is verified (the spherical harmonics form a orthogonal basis for the L 2 
scalar product). 

In the numerical test, we use the P* model for which 
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Remark 5.1. Since the spherical harmonics are eigenvectors of scattering operators of the form, 



Q(f) 



J ^ P (n, n') (/(i, x, si) ~ f(t, x, «)) dri 



or 



Q(/) = A n /(t,x,n), 



(49) 



(50) 



where Afj is the Laplace- Beltrami operator defined on the sphere and p(Q,,Q. ) is an angular repar- 
tition function, therefore the "diffusive - non diffusive " decomposition for P/v models associated to 
the transport equation with these operators is still valid since the assumption (H2) is verified. 



6 Applications to the Sn models 

The Sn models for the transport equation (43) are defined by 
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with fi = /(Hi), tti a discrete direction and u>i and quadrature weight. 

and £> c = | if the velocities are defined in S 2 and D c = | if the velocities are defined in S 1 . Usually 
the quadrature formula is symmetric with respect the rotation of the axis. These systems admit 
the following diffusion limit 

d t E - div (DVE) = (51) 

with E = J\. u>j/j and D = ^ X) j w j^j ® ^j- 

Proposition 6.1. TTie S'jv models can be formulate to 

d t \J + 1 A 1 d x U + -A 2 d y XJ = -~MJ, 

with Uj = *Jwjfj for each j and R = Id — \/w <g) y'w. The vector -y/w is given by the the square 

root of Wj. 

This system satisfies the following properties 

• dimKeri? = 1, 

• A\ and A 2 are diagonals, 

• is an eigenvalue of R with the multiplicity 1 and the eigenvector Ei = (y/wi, ....,y/w n ), 

• I is an eigenvalue of R with the multiplicity n — 1, 

• The matrix R is symmetric with real coefficients, 

• .AiEi = aE 2 , A 2 Ei — aE 3 . 

Proof. We first prove that 1 is an eigenvalue with the multiplicity n — 1. We notice that Id — 
y^w (g) -y/w corresponds to the orthogonal projection on the hyperplane orthogonal to the vector 
y/w. Therefore 1 is the eigenvalue of the matrix R with the multiplicity n — 1. The projection 
in the space generate by -y/w is equal to zero, thus is an eigenvalue of R associated to the 
eigenvector Ei = \/w- I n a second step, we show that the last property of the proposition 6.1 is 
verified. The condition under the quadrature point J^ Wiili = imply that (A\y/w, yw) = an d 
(A 2 y/w, y/w) = 0. Consequently AE 1 G Ker(i?)- L and AE 2 £ Ker(i?)- L . Using 



Eo 



we obtain AiEi = aE 2 , A2E1 = aE^ with a = \/J2i w i^'T = \/J2i w i^T ■ The equality 




J2i w i^i' = \/J2i w i^t comes from 

J~] Wjtoj ® $7, = D c I d , 



with D c = \ or D c = \. U 

Remark 6.2. In dimension one, the Sn models with Gauss-Legendre quadrature, are equivalent 
to the Pn models. 

Remark 6.3. Unlike the case of the Pjy model, the assumptions (H 2 ) are not satisfied for the 
anisotropic scattering (49). 
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The S 2 model used in the numerical examples writes 
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d t v + -A 1 d x u + -A 2 d y xj 
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and the diffusion limit is 



d t E - div 



2cr 7 



with E = ( j ^ ■ Uj). Defining the orthogonal matrix Q and D the diagonal matrix by 
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then in the unknowns V = Q*U the system writes 
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dtV + -A x d x Y 
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7 Numerical results 



(52) 



(53) 



In this section we describe numerical results obtained for the three models, the hyperbolic heat 
equation (equivalent to Pi ) , P 3 and S2 described previously. We give the results for both the diffu- 
sion and the transport regimes. For each model the results are obtained for 3 types of unstructured 
meshes as illustrated in figures 3 and 4. In this section, contour plots are for the first moment of 
the solution that is p = (U, E\) and we recall that E\ is the basis of the kernel of R. 
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Figure 3: Unstructured quadrangular meshes 
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Figure 4: Kershaw mesh 



7.1 The P x model 



We illustrate the behavior of the asymptotic preserving discretization in 2D. We use (24) to solve 
the Pi model and we compare with an exact diffusion solution. 

The test case is based on the fundamental solution of the heat equation with a diffusion coefficient 
equal to one, called SF(t) [EGHOO]. The initial datas are U x (t = 0) = SF(0.01) and U 2 {t = 0)=0, 
U 3 (t = 0)=0. The diffusion solution is given by U x {t) = 5(0.01 + t), U 2 (t) = 0, U 3 (t) = 0. We 
compare the exact diffusion solution, the solution obtained with the scheme without AP corrector, 
the solution obtained with the scheme with AP corrector which admits a non consistent TPFA 
diffusion scheme [BDF11] and the solution obtained with the consistent asymptotic preserving 
scheme (24). The exact diffusion solution is plotted on Cartesian mesh with 150 cells for each 
direction. The numerical solutions are computed on Kershaw mesh with 150 cells for each direction 
and e = 0.0001. When we solve this problem with the classic upwind scheme (fig. 5), we do 





Figure 5: On the left, we plot the first moment p of the exact solution of the test case at the time 
t = 0.01. On the right, we plot p obtained by the classical upwind scheme at the time t = 5 X 10 -5 

not capture correctly the dynamic of the solution. For the TPFA-asymptotic preserving scheme 
(left on fig. 6), the quality of numerical solution is very dependent of the deformation of the mesh 
and the symmetry of the solution is not preserved. The numerical solution given by the nodal AP 
scheme is close to the exact solution. The quality of the numerical solution is not very sensitive to 
the mesh deformations. 
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Figure 6: We plot the first moment p of the solution at t = 0.01 with two different schemes. In the 
left the scheme is a naive multidimensional extension of the usual 1-d AP scheme: this scheme is 
validated only on Delaunay meshes [BDF11]. In the right the scheme is the nodal JL-(b) (24-25) 
asymptotic preserving which is convergent on general meshes. 

7.2 The S 2 model 

We solve the Friedrichs S% system (2) with the matrices (52). The numerical scheme for the 
"diffusive" part is (24)-(25) and we use an implicit time discretization. We define the first moment 

E = Ei w ^ v ^ with u = (0i=i. •••> 0i=4). 

7.2.1 Numerical results in the diffusion regime 

We note SF 2 (t) the fundamental solution of the heat equation with a diffusion coefficient |. At 
the time zero, each unknown Ui is equal at ^^(O.Ol). The solution at the time Tf — 0.01 is 
S*F2(0.01 + Tf). The model is discretized with the JL-(b) nodal scheme for the "diffusive" part 
the upwind scheme for the other part and an implicit time discretization. We obtain the following 
results for the convergence. 
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Table 2: Order of convergence for the S% scheme on Cartesian mesh 
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Table 3: Order of convergence for the S2 scheme on random quadrangular mesh 



The tables (2)-(3)-(5)-(6)-(7) give the order of convergence for some meshes and values of e. In 
the diffusion regime the numerical method converges with the second order. 

These results deserve some remarks. The order of convergence for e — 0.001 and e — 0.0001 de- 
creases a little when the number of cells increase. This phenomena comes from the fact that we 
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Table 4: Order of convergence for the S2 scheme on "smooth" mesh 
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Table 5: Order of convergence for the £2 scheme on Kershaw mesh 
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Table 6: Order of convergence for the S2 scheme on regular triangular mesh 
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Table 7: Order of convergence for the Si scheme on random triangular mesh 



compare the numerical solution of the Si with the exact solution of the diffusion equation. But the 
exact diffusion is an approximation of the 5*2 solution with an error homogeneous to e. Therefore 
when the numerical error is close to e, it is not justified to compare the error numerical with the 
diffusion solution. 



7.2.2 Transport test case 

We verify here that the "diffusive - non diffusive" decomposition and the AP corrector do not 
disturb the convergence in the transport regime. 
Test case 1 

It is a classical transport case. The quantities are initialized by U\ — x\q 4 6] 2 anc ^ Ui — for 
i > 1. We define a = and e = 1. The solution for Ui(t,x) is the initial solution advected with 
the velocity (1,0), the other variables are equal to zero. The final time is Tf = 0.1. Since the 
initial data is discontinuous the theoretical order is 0.5 for the norm L 1 . We show the order for 
the variable E = J^ wJJi in table 8. 

Test case 2 

We note G(x) a Gaussian function. The initial data are given by [7j(x, i = 0) = G(x) and the 
parameters are defined by a = and e = 1. The solution corresponds to the advection of four 
Gaussian functions with advection velocities (0,1), (0,-1), (1,0) et (—1,0). The final time is 0.2. 
We compare the exact and numerical solutions for the quantity E 
the scheme converges with the first order as can seen on figure 7. 



| 5^i=i Ui- For this test case, 
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Table 8: Order of convergence for the S% scheme for the test case 1 





Figure 7: The red curve with square correspond to the numerical error on Cartesian mesh (left) 
and on random mesh (right). The blue curve with circle correspond to the function f . 

Test case 3: 

The initial data is Ui — 6\ t i with 8\\ a Dirac function centered in x = 1 and y = 1. We take 
e = 1 and a = 1. The analytical solution is constructed with 4 Dirac functions advected in each 
direction. We use a random quadrangular mesh. The result is computed using the stabilized-nodal 
scheme (without spurious modes, see [BDF11]) for the "diffusive" part. The result is given by the 
figure (8). 

Remark 7.1. The last test case allows to exhibit a default of the " diffusive-non diffusive" de- 
composition. Indeed the Sn model preserves the positivity of the discrete distribution function 
associated to the linear kinetic equation, consequently all the unknowns are positive. This property 
is not preserve at the discrete level. 

7.3 The P 3 model 

In this subsection we validate our numerical method for the P3 system. We verify the convergence 
in the diffusion limit. After we propose some test cases in the transport regime. 

7.3.1 Numerical results for Pj, in diffusion limit 

Let SFz{t) be the fundamental solution of the heat equation with a diffusion coefficient of |. The 
initial data is U\{t = 0) = Si^O.Ol) and Ui(t = 0) = for i different of zero. The final time is 
Tf — 0.01. The solution, at the final time, is the fundamental solution at t = 0.02. We provide 
convergence order for implicit scheme and semi-implicit scheme obtained using the semi-implicit 
JL-(b) nodal scheme for the "diffusive" part and a modified Rusanov scheme for the other part 
(see subsection 4.2). The time step is given by At — \h 2 with h the step mesh. 
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Figure 8: First moment p of the fundamental solution of 52 model 



Semi-implicit time discretization 
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Implicit time discretization 
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Table 9: Order of convergence for the P3 scheme on Cartesian mesh 



Semi-implicit time discretization 
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Implicit time discretization 
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Table 10: Order of convergence for the P3 scheme on random quadrangular mesh 



The tables (9)-(10)-(12)-(13)-(14) give the order of convergence for some meshes and some 
values of e. The remarks introduced on the convergence results of the asymptotic preserving 
scheme for the 5*2 model are valid for this test case. 
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Semi-Implicit time discretization 
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Implicit time discretization 
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Table 11: Order of convergence for the P 3 scheme on "smooth" mesh 



Semi-implicit time discretization 
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Implicit time discretization 
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Table 12: Order of convergence for the P 3 scheme on Kershaw mesh 



Semi-implicit time discretization 
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Implicit time discretization 
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Table 13: Order of convergence for the P3 scheme on regular triangular mesh 



7.3.2 Fundamental solution for P 3 and P x models 



Now we solve the P3 and P\ systems with XJ\(t = 0) 



>(i,i) 



and Ui(t = 0) = for i ^ 1 [FHSN11]. 



The "diffusive" part is approximated with the JL-(b) nodal scheme. The time discretization is 
implicit. This test case is described in [HMcll]. The final time is T = 1. The exact solution is 
composed of Dirac functions with the velocities A^ (A^ are the eigenvalues of A\n x + A 2 n y ) and 
smooth functions between the Dirac functions. At the beginning the smooth functions are non 
negatives and becomes negatives for large time. For the P\ system, the speed wave is -4= and for 
the P3 system the maximal velocity is approximately 0.86. The numerical results reproduce this 
behavior, see figure 9. 
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Semi-implicit time discretization 
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Implicit time discretization 
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Table 14: Order of convergence for the P3 scheme on random triangular mesh 
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Figure 9: Left the first moment p of the solution fundamental for the P\ model at the time Tf = 1, 
right p for the fundamental solution for the P3 model at the time Tf = 1 

7.3.3 Lattice problem for P 3 and Pi models 

This test case is an example of a complicated geometry. We consider a checker-board with different 
scattering and absorbing opacities on a lattice core (see [SFL11]). It is interesting for neutron 
transport simulations since is a simplification to a reactor core. The geometry is given in the figure 
10. We define a the scattering opacity and a a the absorption opacity. In the black square and 
the striped squares a a — 10 and a = 0. In the white squares a = 1 and a a = 0. The relaxation 
parameter is defined e = 1 in the whole domain. All the unknowns are equal to zero at the time 
0. We solve the Pi and P3 systems with the additional source term 

d t V + 1 Ad x U + -Bd x U = -^PU + S 

e e e z 

where A\, A^ and P are the matrices associated to the Pi or P3 system and Si — —(a a Ui + Q)5n 

with Su a Kronecker product. The source Q = 1 in the black square and Q = in the rest of the 

domain. 

The P3 systems is solved using the JL-(b) scheme for the "diffusive" part. We plot the first moment 

with a logarithmic scale logiQ at the final time Tf = 3.2. The results for Pi and P3 are given for 

the first moment in figures 11. They are the same as those in [SFLll]-[Bru02]. 
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Figure 10: Geometry for the test case, the domain is [0, 7] x [0, 7] 
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Figure 11: In the left, we solve the Pi model and plot the log 10 of the first moment p. In the right, 
we solve the P3 model and plot the log 10 of p. 

8 Conclusion 

We have studied the discretization on distorted meshes of linear hyperbolic systems with stiff source. 
We have proposed a method called "diffusive - non diffusive" decomposition which consists to split 
the hyperbolic system between the hyperbolic heat equation and a other system which is negligible 
in the diffusion regime. Using an asymptotic preserving scheme for the hyperbolic heat equation 
to discretize the "diffusive " part and a classical scheme to discretize the other part, we obtain 
an asymptotic preserving method for the complete system. For the approximation of transport 
equation, we use this decomposition for the simplified models as P/v or Sn approximations. For 
the Pjv systems the decomposition is natural. Since the first and second moments gives the limit 
regime. The others moments are close to e. The high order moments are added only to obtain 
a better approximation in the pure transport regime (er — 0). For the Sjy models, we remark 
that the diagonalized model admits a structure very close to the structure of the Pjv models. 
The "diffusive - non diffusive" decomposition gives consistent schemes for all the regimes. If the 
numerical methods used to discretize the different parts of the decomposition are stable in norm 
L 2 , the method is stable in norm L 2 . Modifying the schemes for the "non diffusive" part we can 
obtain a semi-implicit scheme with a CFL condition independent of e. However this method is 
not optimal for the discretization of Sn models, since our numerical method does not preserve 



2:-; 



the positivity. In the future, it would be interesting to design positive and asymptotic preserving 
method. 
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